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Algorithm for Liapunov Stability Analysis 


An algorithm has been developed for automatic 
computation of a quadratic estimate of the domain 
of stability for the stable equilibrium states of non- 
linear systems of ordinary differential equations. 
Standard linear stability analysis techniques failed to 
predict adequately the stability of the current NASA 
Orbiting Astronomical Observatory coarse-pointing- 
mode control system. 

In developing the algorithm, various minimization 
and random search techniques were utilized to solve 
the min— max problem which yields the estimate; 
gradient search and penalty function techniques were 
found to be generally inapplicable to high order non- 
linear problems. The new methods are the first known 
example of solving a min— max problem via two 
random searches. 

The algorithm is based on Liapunov stability theory, 
the only available sufficiently general approach to 
nonlinear system stability analysis. The primary dif- 
ficulty in application of the theory is the construction 
of an appropriate Liapunov function. This difficulty 
is eliminated by restricting consideration to positive 
definite quadratic-form Liapunov functions. This re- 
striction is substantial, but it results in an estimate 
that is always an ellipsoid in n-space, and is often 
a better estimate than those obtained with more 
complex functions. The fact that the estimate is 
always a hyperellipsoid means it is easier to visualize 
and interpret than other estimates, and some of the 
computations are simplified. 

The algorithm is applicable to systems that are 
described by a set of quasilinear differential equa- 
tions of the form: 

x = Ax + g(x) 

where x is the n vector of state variables, A is a 


stable matrix, g(x) is at least of order x 2 , and x = 0 
is the equilibrium point of interest (i.e., g(0) = 0). 
Most systems that are designed by techniques cur- 
rently used can be treated as indicated above. The 
objective is to determine or estimate the set of initial 
states from which the system will return to the equili- 
brium state x = 0, i.e., to estimate the domain of 
attraction of the equilibrium state. 

The objective is achieved by using LaSalle’s the- 
orem on the extent of asymptotic stability. Basically, 
the theorem states that if a Liapunov function V(x) 
and its total time derivative have certain properties 
within the set of states x such that V(x) < l, 
where l is a constant, then all trajectories beginning 
in tend toward the equilibrium x = 0. The prob- 
lem of finding such a Liapunov function V(x) is 
resolved by restricting consideration to positive defi- 
nite quadratic forms. The largest value of l for which 
the required conditions hold gives the best estimate 
relative to that particular quadratic form. This value 
is found to be the solution of a constrained minimum 
problem. 

In developing the algorithm for application to this 
nonacademic complex problem, it was necessary to 
solve four major computational problems efficiently. 
These problems are: 1) generation of arbitrary posi- 
tive definite matrices; 2) solution of the Liapunov 
matrix equation; 3) solution of a nine- variable con- 
strained minimum problem; 4) solution of a max- 
imization problem involving 45 variables. The first 
was accomplished by developing a parameterization 
of the set of positive definite matrices. The second 
was solved by selecting from the four available tech- 
niques the one with the least error growth with 
increased system dimension. The third was solved 

(continued overleaf) 


This document was prepared under the sponsorship of the National 
Aeronautics and Space Administration. Neither the United States 
Government nor any person acting on behalf of the United States 


Government assumes any liability resulting from the use of the 
information contained in this document, or warrants that such use 
will be free from privately owned rights. 



by devising an efficient random search tailored to 
the geometry of the problem. This was done subse- 
quent to deter minin g that the widely used gradient 
search techniques and penalty function methods are 
totally inapplicable to this problem. The fourth prob- 
lem was partially solved via an “accelerated random 
search.” Again, gradient techniques are totally in- 
applicable to this complex high-dimensional problem. 

The algorithm is shown to be feasible for solution 
of this complex nonlinear stability analysis problem. 
However, the state of the art of search techniques for 
complex high-dimensional problems severely limits 
its immediate application as an analytical or opera- 
tional tool. It is shown that, compared to simulation 
(the only other currently available tool), it promises 
to produce a lower cost solution with specified con- 
fidence or conversely to produce a higher confidence 
result for a given cost 

The algorithm developed and examined in this 
study is a feasible solution to the problem of esti- 
mating the domain of attraction of an equilibrium 
state of a complex nonlinear physical system. Its 
further development into a practical tool is limited, 
however, in two ways. (1) Since the nonlinearities 
are not accounted for in the geometric shape of the 
estimate, it is conceivable that there are systems for 
which one cannot find a suitably large family of 
ellipsoids within the domain of attraction such that 
all trajectories cut them in the inward direction as 
required by the problem formulation. (2) The state 
of the art in search techniques must be substantially 
advanced so that the maximization problem can be 
effectively and efficiently solved. For problems of 
reasonable complexity, gradient techniques must be 
abandoned in favor of random-search techniques. 
The problem is further complicated because the ob- 
jective function is stated in terms of the solution of 
a minimum problem which is dependent upon the 
variables over which the maximization occurs. Thus, 
the problem is one of high-dimension with unknown, 


complex geometry. 

Despite these hurdles, the continued development 
of the algorithm remains attractive because it prom- 
ises to provide a tool that can provide higher con- 
fidence stability information for complex no nlin ear 
systems at a lower cost than conventional simulation 
methods. 
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